module ObjFunction
   use Index
   use OptimisationTypes
   use QSort
   use Parameter
   use Population
   use ObsModel
   use dateutils
   use HTypes
   implicit none

  
contains


recursive function EstimatableParameter(Opt, Cat, iFleet, iYear, iAge) result(Res)
! Purpose: Determine if parameter should be estimated
! In     : Cat - category
!          iFleet - Fleet
!          iYear - Fleet
!          iAge - Fleet
   implicit none
   type(TOptimisation), intent(in) :: Opt
   integer, intent(in) :: Cat, iFleet, iYear, iAge
   logical :: Res
   type(TParamElm), pointer :: pParamElm
   pParamElm => GetParamElm(Opt, Cat, iFleet, iYear, iAge, .false.)
   Res = .false.
   select case(pParamElm%Flag)
      case(flEstimate) 
         Res = .true.
      case(flSpecialCode) 
         Res = Opt%IndexInput%ModelType /= pmtVPA
         ! Separable and ISVPA models threat special code as 1.
         ! Only vpa models use special code for internal F calculation.
!      case (flEqualPrevYear)
!         if (0 < iYear) then
!            Res = EstimatableParameter(Opt, Cat, iFleet, iYear - 1, iAge)
!         end if
!      case (flEqualPrevAge)
!         if (0 < iAge) then
!            Res = EstimatableParameter(Opt, Cat, iFleet, iYear, iAge - 1)
!         end if
   end select 
end function EstimatableParameter


subroutine SearchForActiveParameters(Opt, DoAssociate, N)
! Purpose : Search for active parameters
!           Active parameters are estimated based upon 2 conditions:
!           1. Relevance according to model type 
!           2. Estimate flag (=1) 
! In      : DoAssociate - true if pactive should be initialized
! Out     : N - number of active parameters
!   History: When        Who        What
!            07-06-07    Skalevik   Created
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   logical, intent(in) :: DoAssociate
   integer, intent(out) :: N
   integer :: Cat, iFleet, iYear, iAge, NFleets, NYears, NAges, PrevFlag
   type(TParamCube), pointer :: pParamCube
   type(TParamElm), pointer :: pParamElm
   logical :: Estimatable, Optimisable
   ! Calculate the number of parameters to estimate
   N = 0
   PrevFlag = -1
   do Cat = cFirstCategory, cLastCategory
      pParamCube => Opt%ParamCubes(Cat)
      NFleets = size(pParamCube%Elm, 1)
      do iFleet = 0, NFleets - 1
         NYears = min(NumYears(Opt%IndexInput), size(pParamCube%Elm, 2))
         !NYears = size(pParamCube%Elm, 2)
         do iYear = 0, NYears - 1
            NAges = size(pParamCube%Elm, 3)
            do iAge = 0, NAges - 1
               pParamElm => GetParamElm(Opt, Cat, iFleet, iYear, iAge, .false.)
               Estimatable = .true.
               ! Check if parameter is Estimatable
               select case (Opt%IndexInput%ModelType)
                  case (pmtSep)
                     select case (Cat)
                        case (cNValues)
                           Estimatable = (iYear == NYears - 1) .or. (iAge == NumTrueAges(Opt%IndexInput) - 1)
                     end select
                  case (pmtVPA)
                     select case (Cat)
                        case (cNValues)
                           !Estimatable = (iYear == NumCatchYears(Opt%IndexInput) - 1) .and. (iAge == NumTrueAges(Opt%IndexInput) - 1)
                           ! N values are estimatable in the terminal projection area at youngest age
                           Estimatable = (iYear == NumCatchYears(Opt%IndexInput) - 1) .or. &
                              ((iAge == NumTrueAges(Opt%IndexInput) - 1) .and. (iYear <= NumCatchYears(Opt%IndexInput) - 1)) .or. &
                              ((iAge == 0) .and. (iYear > NumCatchYears(Opt%IndexInput) - 1))
                        case(cSelection, cFYearFactor) 
                           Estimatable = .false.
                        !case default
                     end select
                  case (pmtISVPA)
                     select case (Cat)
                        case (cNValues)
                           ! N values are estimatable in the terminal projection area at youngest age
                           Estimatable = (iAge == 0) .and. (iYear > NumCatchYears(Opt%IndexInput) - 1)
                        case (cSelection)
                           Estimatable = .false.
                        case (cFYearFactor)
                           Estimatable = iYear == NumCatchYears(Opt%IndexInput) - 1
                     end select
               end select
               
               ! 2 Resolve the estimated flag 
               Optimisable = Estimatable .and. EstimatableParameter(Opt, Cat, iFleet, iYear, iAge) 
               ! Assosiate a pactive element (pointer to double precision) to a TParamElm%Value
               
               if (Optimisable) then
                  if (DoAssociate) then
                     Opt%pactive(N) = pParamElm
                  end if
                  N = N + 1
               end if
               select case(pParamElm%Flag)
                  case(flKeepFixed, flEstimate) 
                     PrevFlag = pParamElm%Flag
               end select
            end do
         end do
      end do
   end do
end subroutine SearchForActiveParameters


subroutine InitializeDerivedParameters(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: Cat, iFleet, iYear, iAge, NFleets, NYears, NAges
   type(TParamCube), pointer :: pParamCube
   type(TParamElm), pointer :: pParamElm
   ! Initialize all derived parameters 
   do Cat = cFirstCategory, cLastCategory
      pParamCube => Opt%ParamCubes(Cat)
      NFleets = size(pParamCube%Elm, 1)
      do iFleet = 0, NFleets - 1
         NYears = size(pParamCube%Elm, 2)
         do iYear = 0, NYears - 1
            NAges = size(pParamCube%Elm, 3)
            do iAge = 0, NAges - 1
               pParamElm => GetParamElm(Opt, Cat, iFleet, iYear, iAge, .false.)
               select case(pParamElm%Flag)
                  case (flKeepFixed, flEstimate)
                     call SetPValueEx(pParamElm, pParamElm%value)
               end select
            end do
         end do
      end do
   end do
end subroutine InitializeDerivedParameters


subroutine InitializeActiveParameters(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: N
   call SearchForActiveParameters(Opt, .false., N)
   if (N == 0) return
   allocate(Opt%pactive(0 : N - 1))
   call SearchForActiveParameters(Opt, .true., N)
   call InitializeDerivedParameters(Opt)
end subroutine InitializeActiveParameters


!function GetNumPhiValues(Opt)
!   implicit none
!   type(TOptimisation) :: Opt
!   integer :: GetNumPhiValues, iDataSourceType, DataSourceType, iDataSource, NumDataSources
!   type(TObservedData), pointer :: pObsData
!   GetNumPhiValues = 0
!   do iDataSourceType = 0, size(c_DataSource) - 1
!      DataSourceType = c_DataSource(iDataSourceType)
!      NumDataSources = GetNumObservedDataElms(Opt, DataSourceType)
!      if (0 == NumDataSources) cycle
!      do iDataSource = 0, NumDataSources - 1
!         pObsData => Opt%ObservedData(DataSourceType)%Elm(iDataSource)
!         GetNumPhiValues = GetNumPhiValues + pObsData%NumYears * pObsData%NumAges
!      enddo
!   enddo
!end function GetNumPhiValues


subroutine FillInM(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: iAge, iYear, NYears, NAges
   NYears = NumYears(Opt%IndexInput) + 1
   NAges = NumAges(Opt%IndexInput)
   do iYear = 0, NYears - 1
      do iAge = 0, NAges - 1
         if (iYear == NYears - 1) then ! intermediate year
            Opt%PopModel%May(iYear, iAge) = 0.0d0
         else
            Opt%PopModel%May(iYear, iAge) = PValue(Opt, cNatMort, 0, iYear, iAge) 
         endif
      enddo
   enddo
end subroutine FillInM


subroutine ResetOpt(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   double precision :: o
   o = objfunc(Opt, .false.)
   Opt%Loc%ObjFunc = o
   Opt%Loc%Iter = 0
   Opt%Loc%Delta = 0.25
   Opt%Loc%Improved = .true.
   Opt%seed = 15349668.0
end subroutine ResetOpt

   
function objfunc(Opt, Finished)
! Purpose : Calculates the object function value, the optimisation fitness
!
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   double precision :: objfunc
   integer :: iDataSourceType, DataSourceType, iDataSource, iYear, iAge, iObsYear, iObsAge, I, &
      IStart, iFrom, iTo, NumDataSources, NYears, M
   type(TDataSource), pointer :: pDataSource
   type(TObservedData), pointer :: pObsData
   double precision :: ModVal, ObsVal, Phi, PhiWeight, Res, T, objval, k, My, Med
   logical :: Finished
   logical :: IgnorePhi
   !type(TObservedData) :: Obs
   !open(90, file = GetFileName(Opt%ProjectDir, "ObjFunc"), action = 'write', status =  'replace')
   objfunc = 0
   NYears = Opt%IndexInput%LastYear - Opt%IndexInput%FirstYear + 1
   if ((0 == NYears) .or. (0 == NumAges(Opt%IndexInput))) then
      return
   end if
   call CalculatePopulation(Opt, Finished)
   if(Finished) then
      if (Opt%IndexInput%ModelType == pmtISVPA) then
         ! Scale the selection and year factor
         !call NormalizeSelection(Opt, .true.)
      endif
      call FillInM(Opt)
   endif

!   if (Finished) then
!      open(cPhiFile, file = GetWorkFileName(Opt%ProjectDir, s_PhiFile), action = 'write', status = 'replace')
!      write(cPhiFile,*) ' Idx    Yr    Age     Obs       Mod       Wgt       Phi   '   
!   end if
   ! Loop: For each optimisable data source 
   do iDataSourceType = 0, size(c_DataSource) - 1
      DataSourceType = c_DataSource(iDataSourceType)
      NumDataSources = GetNumObservedDataElms(Opt, DataSourceType)
      ! Do not handle if data source are not used
      if (0 == NumDataSources) cycle

      ! VPA doesnt use catch residuals
      !IgnorePhi = (DataSourceType == odCanum) .and. & 
      !((Opt%IndexInput%ModelType == pmtVPA) .or. (Opt%IndexInput%ModelType == pmtISVPA))

      ! Loop: For each index in data source
      do iDataSource = 0, NumDataSources - 1
         pDataSource => GetDataSource(Opt, DataSourceType, iDataSource)
         pObsData => Opt%ObservedData(DataSourceType)%Elm(iDataSource)
         I = 0
         IStart = 0
         !Obs = GetObsData(Opt, DataSourceType, iDataSource)
         ! Loop: For each year 
         
         !pObsData%PhiArray(0 : pObsData%NumYears * pObsData%NumAges - 1) = 0.0d0
         !if (Finished) then
         !   pObsData%Phi(0 : pObsData%NumYears - 1, 0 : pObsData%NumAges - 1) = 0.0d0
         !   pObsData%Res(0 : pObsData%NumYears - 1, 0 : pObsData%NumAges - 1) = 0.0d0
         !   pObsData%Cbl(0 : pObsData%NumYears - 1, 0 : pObsData%NumAges - 1) = 0.0d0
         !   pObsData%ModVal(0 : pObsData%NumYears - 1, 0 : pObsData%NumAges - 1) = 0.0d0
         !endif
         pObsData%PhiArray(0 : size(pObsData%PhiArray, dim = 1) - 1) = 0.0d0
         pObsData%ResArray(0 : size(pObsData%ResArray, dim = 1) - 1) = 0.0d0
         pObsData%NumPositiveUnweightedResiduals = 0
         if (Finished) then
            pObsData%Phi(0 : size(pObsData%Phi, dim = 1) - 1, 0 : size(pObsData%Phi, dim = 2) - 1) = 0.0d0
            pObsData%ResU(0 : size(pObsData%ResU, dim = 1) - 1, 0 : size(pObsData%ResU, dim = 2) - 1) = 0.0d0
            pObsData%ResW(0 : size(pObsData%ResW, dim = 1) - 1, 0 : size(pObsData%ResW, dim = 2) - 1) = 0.0d0
            pObsData%CblFac(0 : size(pObsData%CblFac, dim = 1) - 1, 0 : size(pObsData%CblFac, dim = 2) - 1) = 0.0d0
            pObsData%CblPow(0 : size(pObsData%CblPow, dim = 1) - 1, 0 : size(pObsData%CblPow, dim = 2) - 1) = 0.0d0
            pObsData%NFromDS(0 : size(pObsData%NFromDS, dim = 1) - 1, 0 : size(pObsData%NFromDS, dim = 2) - 1) = 0.0d0
            pObsData%ModVal(0 : size(pObsData%ModVal, dim = 1) - 1, 0 : size(pObsData%ModVal, dim = 2) - 1) = 0.0d0
         endif

         !if(pDataSource%Weight == 0.0d0) cycle
         IgnorePhi = pDataSource%Weight == 0.0d0

         do iObsYear = 0, pObsData%NumYears - 1
            do iObsAge = 0, pObsData%NumAges - 1
               iYear = iObsYear + pObsData%FirstYear - Opt%IndexInput%FirstYear
               iAge = iObsAge + pObsData%FirstAge - Opt%IndexInput%FirstAge 
               if ((iYear < 0) .or. (iYear > NYears - 1) .or. &
		            (iAge < 0) .or. (iAge > NumAges(Opt%IndexInput) - 1)) then
		            cycle
		         endif
               ! Action GetObsDataValue: Retrieve observed data value
               ObsVal = GetObsDataValue(Opt, DataSourceType, iDataSource, iYear, iAge)
               ! Reminder: Is Phi set to 0 if Obs = 0???
               ! If: missing observation   
               Phi = 0
               ModVal = 0.0d0
               PhiWeight = GetObsDataValue(Opt, DataSourceType + 1, iDataSource, iYear, iAge) 
               PhiWeight = PhiWeight * pDataSource%Weight
               if (ObsVal >= 0) then
                  ! Action GetObsModDataValue: Retrieve modeled data value 
                  ModVal = GetObsModDataValue(Opt, DataSourceType, iDataSource, iYear, iAge, .true.)
                  if ((.not. IgnorePhi) .and. (ModVal >= 0)) then
                     select case(pDataSource%ObjType)
                        case(otSSQ, otMedian, otAMD)
                           My = pDataSource%My
                           if (((ObsVal + My) > 0) .and. ((ModVal + My) > 0)) then
                              Phi = log((ObsVal + My) / (ModVal + My)) ** 2
                           endif
                        case(otGammaLikelihood)
                           if ((ObsVal > 0) .and. (ModVal > 0)) then
                              k = PValue(Opt, GammaCategoryFromDataSource(DataSourceType), iDataSource, 0, 0)
                              Phi = -((k - 1) * log(ObsVal) - k * ObsVal / ModVal - k * log(ModVal / k) - StirlingApprox(k))
                           endif
                        case(otGammaDeviance)
                           if ((ObsVal > 0) .and. (ModVal > 0)) then
                              Phi = -log(ObsVal / ModVal) + (ObsVal - ModVal) / ModVal
                           endif
                        case(otPoisson)
                           if ((ObsVal > 0) .and. (ModVal > 0)) then
                              Phi = ModVal - ObsVal + ObsVal * log(ObsVal / ModVal)
                           else
                              Phi = ModVal 
                           endif
                     end select
                     ! Action ObjFunc_SingleValue: Retrieve phi value
                     if (0 <= PhiWeight) then
                        Phi = Phi * PhiWeight
                     end if
                  endif
               endif
               if (Finished) then
                  if ((ObsVal > 0) .and. (ModVal > 0)) then
                     Res = log(ObsVal / ModVal)
                  else
                     Res = 0
                  end if
                  call SetResDataValue(Opt, DataSourceType, iDataSource, iYear, iAge, ModVal, Phi, PhiWeight, Res);
               end if
               pObsData%PhiArray(I) = Phi 
               I = I + 1
            end do
         end do    
         
         if(IgnorePhi) cycle
         select case(pDataSource%ObjType)
            case(otMedian, otAMD) 
               ! Median takes M values from the middle of a sorted array
               call QuickSort(pObsData%PhiArray, 0, size(pObsData%PhiArray) - 1, size(pObsData%PhiArray))

               IStart = 0
               do IStart = 0, size(pObsData%PhiArray) - 1
                  if (pObsData%PhiArray(IStart) > 0) exit
               enddo
               M = Opt%IndexInput%NumMedianElms
               T = IStart + (size(pObsData%PhiArray) - IStart) / 2
               iFrom = max(IStart, int(anint(T - M / 2)))
               iTo = min(size(pObsData%PhiArray) - 1, iFrom + M - 1)
               
               if (pDataSource%ObjType == otAMD) then
                  Med = 0
                  do i = iFrom, iTo
                     Med = Med + pObsData%PhiArray(i)
                  enddo
                  Med = Med / (iFrom - iTo + 1)
                  do i = 0, size(pObsData%PhiArray) - 1
                     if (pObsData%PhiArray(i) /= 0) then
                        pObsData%PhiArray(i) = abs(pObsData%PhiArray(i) - Med)
                     endif
                  enddo
                  call QuickSort(pObsData%PhiArray, 0, size(pObsData%PhiArray) - 1, size(pObsData%PhiArray))
                  IStart = 0
                  do IStart = 0, size(pObsData%PhiArray) - 1
                     if (pObsData%PhiArray(IStart) > 0) exit
                  enddo
                  M = Opt%IndexInput%NumMedianElms
                  T = IStart + (size(pObsData%PhiArray) - IStart) / 2
                  iFrom = max(IStart, int(anint(T - M / 2)))
                  iTo = min(size(pObsData%PhiArray) - 1, iFrom + M - 1)
               endif
            case default
               iFrom = 0
               iTo = size(pObsData%PhiArray) - 1
         end select
         objval = 0
         do i = iFrom, iTo
            objval = objval + pObsData%PhiArray(i)
         enddo
         objfunc = objfunc + objval
      end do
   end do
end function objfunc


subroutine InitializeRun(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   ! Reallocate run
   call DeallocateRun(Opt)
   call InitializeActiveParameters(Opt)
   call Population_Allocate(Opt)
   allocate(Opt%iseq(size(Opt%pactive)))
   allocate(Opt%seq(size(Opt%pactive)))
   ! After the run is initialized, reset the run:
   call ResetOpt(Opt)
end subroutine InitializeRun


subroutine DeallocateRun(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   if (associated(Opt%iseq)) deallocate(Opt%iseq)
   if (associated(Opt%seq)) deallocate(Opt%seq)
   if (associated(Opt%pactive)) deallocate(Opt%pactive)
   call DeallocatePopulation(Opt)
end subroutine DeallocateRun


subroutine CreateNewProject(Opt, FirstYear, LastYear, FirstAge, LastAge, NumFleets, NumSSB, NumBiomass)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: FirstYear, LastYear, FirstAge, LastAge, NumFleets, NumSSB, NumBiomass
   ! Reallocate run
   call CloseRun(Opt)
   ! Read index, data and parameters from file:
   allocate(Opt%Loc)
   call CreateNewIndexData(Opt, FirstYear, LastYear, FirstAge, LastAge, NumFleets, NumSSB, NumBiomass)
   call CreateNewObsData(Opt)
   ! Reallocate parameter structure
   call DeallocateParameters(Opt)
   call InitializeParameterStructures(Opt)
   call InitializeRun(Opt)
   Opt%Loc%Stage = 0
end subroutine CreateNewProject

   
function LoadRun(Opt)
   implicit none
   logical :: LoadRun
   type(TOptimisation), intent(inout) :: Opt

   Opt%ErrorStr = ""
   Opt%RunErrorStr = ""
   LoadRun = .false.
   ! Reallocate run
   call CloseRun(Opt)
   ! Read index, data and parameters from file:
   allocate(Opt%Loc)
   if(.not. ReadIndexData(Opt)) return
   if(.not. ReadObsData(Opt)) return
   if(.not. LoadParameters(Opt)) return
   ! Initialize the run structure:
   if(Opt%RunErrorStr == "") then
      call InitializeRun(Opt)
   endif

   ! After run is loaded and intialized, do following:
   ! Reset the stage
   Opt%Loc%Stage = 0
   !Reset the run id
   Opt%RunId = datetimestr()
   ! Clear message log file
   open(cMessageFile, file = GetWorkFileName(Opt, s_LogFile), status = 'replace')
   close(cMessageFile)
   ! Clear debug tuning file
   if(Opt%DumpTuningDetails) then
      open(cTuningFile, file = GetWorkFileName(Opt, s_TuningFile), status = 'replace')
         write(cTuningFile,*) ' Idx    Delta   Value'   
      close(cTuningFile)
   endif
   LoadRun = .true.
end function LoadRun


subroutine CloseRun(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   
   call DeallocateRun(Opt)
   call DeallocateParameters(Opt)
   call DeallocateObsData(Opt)
   call DeallocateIndexData(Opt)
   if (associated(Opt%Loc)) deallocate(Opt%Loc)
end subroutine CloseRun


subroutine CalculateStockSummary(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   double precision :: FRef, SSB, Biomass
   integer :: iYear
   double precision :: obj

   obj = objfunc(Opt, .true.)

   do iYear = 0, NumYears(Opt%IndexInput)  
      Opt%PopModel%StockSummary(iYear, ssRecruits) = Opt%PopModel%Nay(iYear, 0)
      SSB = GetObsModDataValue(Opt, odSSB, 0, iYear, 0, .false.)
      Opt%PopModel%StockSummary(iYear, ssSSB) = SSB
      Biomass = GetObsModDataValue(Opt, odBiomass, 0, iYear, 0, .false.)
      Opt%PopModel%StockSummary(iYear, ssTotalBiomass) = Biomass
      Opt%PopModel%StockSummary(iYear, ssFMean) = GetFRef(Opt, iYear)
      Opt%PopModel%StockSummary(iYear, ssLandings) = GetObsDataValue(Opt, odCaTon, 0, iYear, 0)
   enddo
   
end subroutine CalculateStockSummary


subroutine SavePopulationModel(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: ioerr, DataSourceType, iDataSource, iDataSourceType, NumDataSources, iYear, iFirstYear, iLastYear 
   character(len = 256) :: S 
   character(len = 256) :: FileName
   double precision :: r
   integer :: NumDig, NumDec
   type(TParamElm), pointer :: pParamElm
   type(TObservedData) :: Obs
   call CalculateStockSummary(Opt)

   FileName = Opt%IndexInput%OutFiles(roSummary)
   if ((Opt%IndexInput%ProcessType == ptRetrospectiveRun) .and. &
      (Opt%Loc%Stage > 0)) then
      FileName = "Summary.retro" // lrtrim(inttostr(Opt%Loc%Stage)) // ".txt"
   endif
   FileName = GetWorkFileName(Opt, lrtrim(FileName))
   open(cFileUnit, file = FileName, action = 'write', iostat = ioerr, &
      status = 'replace')
   if (io_no_error .ne. ioerr) then
      return
   endif
   write(cFileUnit, "(A)") "Summary output"
   write(cFileUnit, "(A)") "Run id: " // lrtrim(Opt%RunId)
   write(cFileUnit, "(A)") "Process: " // lrtrim(ProcessTypeDescr(Opt%IndexInput%ProcessType))
   write(cFileUnit, "(A)") "Model: " // lrtrim(ModelTypeDescr(Opt%IndexInput%ModelType, &
      Opt%IndexInput%ISVPAControlType))
   write(cFileUnit,"(A)") ""
   write(cFileUnit, '(A)') 'N in model:'
   call print_matrix(Opt%PopModel%Nay, out_unit=cFileUnit, col_start = Opt%IndexInput%FirstYear, &
      row_start = Opt%IndexInput%FirstAge, &
      col_width = 18, &
      num_dec = 4 &
      )
   write(cFileUnit, *)

   write(cFileUnit, '(A)') 'F in model:'
   call print_matrix(Opt%PopModel%Fay, out_unit=cFileUnit, col_start = Opt%IndexInput%FirstYear, &
      row_start = Opt%IndexInput%FirstAge, &
      col_width = 18, &
      num_dec = 4 &
      )
   write(cFileUnit, *)

   write(cFileUnit, '(A)') 'M in model:'
   call print_matrix(Opt%PopModel%May, out_unit=cFileUnit, col_start = Opt%IndexInput%FirstYear, &
      row_start = Opt%IndexInput%FirstAge, &
      col_width = 18, &
      num_dec = 4 &
      )
   write(cFileUnit, *)

   write(cFileUnit, '(A)') 'Summary:'
   write(S, '(A1I2A1I2)') "F", Opt%IndexInput%FirstRefAge, "-", Opt%IndexInput%LastRefAge
   write(cFileUnit, '(A4, A15, A15, A15, A15, A15)') 'Year', 'Recruit', 'TotBiomass', 'SSB', lrtrim(S), 'Landings'
   do iYear = 0, NumYears(Opt%IndexInput)  
      !write(cFileUnit, '(I4, F15.4, F15.3, F15.3, F15.3, F15.3)') &
      !Opt%IndexInput%FirstYear + iYear, &
      !Opt%PopModel%StockSummary(iYear, ssRecruits), &
      !Opt%PopModel%StockSummary(iYear, ssTotalBiomass), &
      !Opt%PopModel%StockSummary(iYear, ssSSB), &
      !Opt%PopModel%StockSummary(iYear, ssFMean), &
      !Opt%PopModel%StockSummary(iYear, ssLandings)
      write(cFileUnit, '(I4)', advance='no') Opt%IndexInput%FirstYear + iYear
      call WriteFToFile(cFileUnit, Opt%PopModel%StockSummary(iYear, ssRecruits), 15, 4, .false.)
      call WriteFToFile(cFileUnit, Opt%PopModel%StockSummary(iYear, ssTotalBiomass), 15, 3, .false.)
      call WriteFToFile(cFileUnit, Opt%PopModel%StockSummary(iYear, ssSSB), 15, 3, .false.)
      call WriteFToFile(cFileUnit, Opt%PopModel%StockSummary(iYear, ssFMean), 15, 3, .false.)
      call WriteFToFile(cFileUnit, Opt%PopModel%StockSummary(iYear, ssLandings), 15, 3, .true.)
   enddo
   close(cFileUnit);

   FileName = Opt%IndexInput%OutFiles(roOutput)
   if ((Opt%IndexInput%ProcessType == ptRetrospectiveRun) .and. &
      (Opt%Loc%Stage > 0)) then
      FileName = "Output.retro" // lrtrim(inttostr(Opt%Loc%Stage)) // ".txt"
   endif
   FileName = GetWorkFileName(Opt, lrtrim(FileName))
   open(cDsFile, file = FileName, action = 'write', iostat = ioerr, &
      status = 'replace')
   write(cDsFile, "(A)") s_CurrentVersion
   write(cDsFile, "(A)") "Data source output"
   write(cDsFile, "(A)") "Run id: " // lrtrim(Opt%RunId)
   write(cDsFile, "(A)") ""
   do iDataSourceType = 0, size(c_DataSource) - 1
      DataSourceType = c_DataSource(iDataSourceType)
      NumDataSources = GetNumObservedDataElms(Opt, DataSourceType)
      if (0 == NumDataSources) cycle
      do iDataSource = 0, NumDataSources - 1
         Obs = GetObsData(Opt, DataSourceType, iDataSource)
         
         write(cDsFile, '(A)') GetDataSourceTechName(DataSourceType, iDataSource)
         write(cDsFile, '(A)') Obs%Title 

         write(cDsFile, '(A)') 'Observed values'
         iFirstYear = max(0, Opt%IndexInput%FirstYear - Obs%FirstYear)
         iLastYear = min(Obs%NumYears - 1, Opt%IndexInput%LastYear - Obs%FirstYear)
         call print_matrix(Obs%Obs(iFirstYear : iLastYear, :), out_unit=cDsFile, col_start = &
            Obs%FirstYear + iFirstYear, &
            row_start = Obs%FirstAge, &
            col_width = GetOutFormatColWidth(Obs%NumDig, Obs%NumDec + 2), &
            num_dec = Obs%NumDec + 2)
            !fmt_row = "(i6,1000" // lrtrim(GetOutFormat(Obs%NumDig, Obs%NumDec + 2)) // ")", &
            !fmt_col = "(6x,1000i" // lrtrim(inttostr(3 + Obs%NumDig + 1 + Obs%NumDec + 2)) // ")")
         write(cDsFile, '(A)')

         write(cDsFile, '(A)') 'Modelled values'
         call print_matrix(Obs%ModVal(iFirstYear : iLastYear, :), out_unit=cDsFile, col_start = &
            Obs%FirstYear + iFirstYear, &
            row_start = Obs%FirstAge, &
            col_width = GetOutFormatColWidth(Obs%NumDig, Obs%NumDec + 2), &
            num_dec = Obs%NumDec + 2)
            !fmt_row = "(i6,1000" // lrtrim(GetOutFormat(Obs%NumDig, Obs%NumDec + 2)) // ")", &
            !fmt_col = "(6x,1000i" // lrtrim(inttostr(3 + Obs%NumDig + 1 + Obs%NumDec + 2)) // ")")
         write(cDsFile, '(A)')
         
         write(cDsFile, '(A)') 'Phi values'
         call print_matrix(Obs%Phi(iFirstYear : iLastYear, :), out_unit=cDsFile, col_start = &
            Obs%FirstYear + iFirstYear, &
            row_start = Obs%FirstAge, &
            col_width = GetOutFormatColWidth(Obs%NumDig, Obs%NumDec + 5), &
            num_dec = Obs%NumDec + 5, &
            !fmt_row = "(i6,1000" // lrtrim(GetOutFormat(Obs%NumDig, Obs%NumDec + 5)) // ")", &
            !fmt_col = "(6x,1000i" // lrtrim(inttostr(3 + Obs%NumDig + 1 + Obs%NumDec + 5)) // ")", &
            summarize = .true.)
         write(cDsFile, '(A)')
         
         if (DataSourceType == odCanum) then
            write(cDsFile, '(A)') 'Selection'
            ! Canum selection stored in CblFac
            call print_matrix(Obs%CblFac(iFirstYear : iLastYear, :), out_unit=cDsFile, col_start = &
               Obs%FirstYear + iFirstYear, &
               row_start = Obs%FirstAge, &
               col_width = 16, &
               num_dec = 6 &
               )
            ! Canum does not have catchability power!
         else
            write(cDsFile, '(A)') 'Catchability factors'
            call print_matrix(Obs%CblFac(iFirstYear : iLastYear, :), out_unit=cDsFile, col_start = &
               Obs%FirstYear + iFirstYear, &
               row_start = Obs%FirstAge, &
               col_width = GetOutFormatColWidth(Obs%NumDig, Obs%NumDec + 2), &
               num_dec = Obs%NumDec + 2)
            write(cDsFile, '(A)')
            write(cDsFile, '(A)') 'Catchability powers'
            call print_matrix(Obs%CblPow(iFirstYear : iLastYear, :), out_unit=cDsFile, col_start = &
               Obs%FirstYear + iFirstYear, &
               row_start = Obs%FirstAge, &
               col_width = GetOutFormatColWidth(Obs%NumDig, Obs%NumDec + 2), &
               num_dec = Obs%NumDec + 2)
         end if
         write(cDsFile, '(A)')
         

         write(cDsFile, '(A)') 'Weightings'
         !Obs = GetObsData(Opt, DataSourceType + 1, iDataSource)
         call print_matrix(Obs%PhiWeight(iFirstYear : iLastYear, :), out_unit=cDsFile, col_start = &
            Obs%FirstYear + iFirstYear, &
            row_start = Obs%FirstAge, &
            col_width = GetOutFormatColWidth(Obs%NumDig, Obs%NumDec + 5), &
            num_dec = Obs%NumDec + 5)
            !fmt_row = "(i6,1000" // lrtrim(GetOutFormat(Obs%NumDig, Obs%NumDec + 5)) // ")", &
            !fmt_col = "(6x,1000i" // lrtrim(inttostr(3 + Obs%NumDig + 1 + Obs%NumDec + 5)) // ")")
         write(cDsFile, '(A)')

         write(cDsFile, '(A)') 'Unweighted residuals'
         !Obs = GetObsData(Opt, DataSourceType + 1, iDataSource)
         call print_matrix(Obs%ResU(iFirstYear : iLastYear, :), out_unit=cDsFile, col_start = &
            Obs%FirstYear + iFirstYear, &
            row_start = Obs%FirstAge, &
            col_width = GetOutFormatColWidth(10, 4), &
            num_dec = 4)
         write(cDsFile, '(A)')

         write(cDsFile, '(A)') 'Weighted residuals'
         !Obs = GetObsData(Opt, DataSourceType + 1, iDataSource)
         call print_matrix(Obs%ResW(iFirstYear : iLastYear, :), out_unit=cDsFile, col_start = &
            Obs%FirstYear + iFirstYear, &
            row_start = Obs%FirstAge, &
            col_width = GetOutFormatColWidth(10, 4), &
            num_dec = 4)
         write(cDsFile, '(A)')

         write(cDsFile, '(A)') 'N from datasource'
         !Obs = GetObsData(Opt, DataSourceType + 1, iDataSource)
         call print_matrix(Obs%NFromDS(iFirstYear : iLastYear, :), out_unit=cDsFile, col_start = &
            Obs%FirstYear + iFirstYear, &
            row_start = Obs%FirstAge, &
            col_width = GetOutFormatColWidth(10, 4), &
            num_dec = 4)
         write(cDsFile, '(A)')
      end do
   end do
   close(cDsFile)
end subroutine SavePopulationModel



end module ObjFunction
